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SUMMARY 


A three-dimensional laminar and turbulent boundary-layer method is devel- 
oped for compressible flow over swept wings. The governing equations and cur- 
vature terms are derived in detail for a nonorthogonal, curvilinear coor 1 - 
nate system. Reynolds shear-stress terms are modeled by the Cebeci-Smith 
eddy-viscosity formulation. The governing equations are descntized using the 
second-order accurate, predictor-corrector finite-difference technique of 
Matsuno, which has the advantage that the crossflow difference formulas are 
formed independent of the sign of the crossflow velocity component. 


The method is coupled with a full potential wing-body inviscid code 
(FLO-30) and the inviscid-viscous interaction is performed by updating the 
oriqinal winq surface with the viscous displacement surface calculated by the 
boundary- layer code. The number of these "global" iterations ranged from five 
to twelve depending on the Mach number, sweep angle, and angle of attack. 
Several test cases are computed by this method and the results are compared 
with another inviscid-viscous interaction method (TAWFIVE) and with 
experimental data. 
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INTRODUCTION 


To obtain accurate flow field calculations over wings and wing-body 
combinations, the effects of viscosity must be accounted for in the 
solution. This may be accomplished by either solving some form of the Navier- 
Stokes equations or by interacting an inviscid code with a viscous boundary- 
layer method. The solution of the three-dimensional Navier-Stokes equations 
requires extremely large computer storage in order to resolve the various 
length scales and thus, also requires very long calculation run times. 
Therefore a great deal of attention has been focused on the more economical 
inviscid-viscous interaction procedure. 

The interaction is performed by updating the inviscid solution by either 
a displacement surface or transpiration velocity component calculated by the 
boundary- layer code and then rerunning the inviscid code. This procedure is 
repeated until successive viscous updates are less than some specified 
tolerance. The methods of Streett [1] and McLean and Randall [2] interact a 
full potential inviscid code with a three-dimensional boundary- layer method by 
use of the displacement surface. A code in which the interaction is performed 
by the transpiration velocity component was developed by Wai, et al. [3]. 

Each of these methods requires the solution of some form of the three- 
dimensional boundary-layer equations. 

Streett [1] and Wai [3] chose to solve the boundary- layer equations in 
the integral form, where the governing partial differential equations have 
been integrated in the normal direction. Swafford and Whitfield [4] developed 
a time-dependent integral method for the three-dimensional boundary- layer 
equations but did not perform the inviscid-viscous interaction. A drawback of 
the integral approach is that one must choose, a priori, the shape of the 
velocity profiles. For example, in Streett' s method the crossflow velocity 
profile cannot predict a crossflow profile with velocity crossover, a 
situation which occurs quite often for three-dimensional flows. Also, 
integral methods are generally less flexible than differential methods. 

The method of McLean and Randall [2] solves the boundary- layer equations 
in the differential form for a body-oriented curvilinear coordinate system. 
However, the coordinate system they used is orthogonal with one surface 
coordinate along the constant chord lines and the other perpendicular to 
them. This type of coordinate system requires a large clustering of grid 
points in the trailing edge region to resolve the outboard portion of a swept 
and tapered wing. Further, the wing root- and tip-chord lines do not lie 
along grid lines so the solution may not be determined along these lines. The 
use of a nonorthogonal coordinate system alleviates these problems. 

Viscous methods which solve the three-dimensional boundary-layer equa- 
tions in a nonorthogonal coordinate system have been developed by Cebeci, 
et al. [5], Matsuno [6], and Van Dalsem and Steger [7]. However, none of 
these methods were interacted with an inviscid code. In this paper, a three- 
dimensional laminar and turbulent boundary- layer method is developed for 
compressible flow over arbitrary swept wings and coupled to a full potential 
inviscid code through a displacement surface calculation. The inviscid 
solution procedure used to drive the boundary-layer method is the FLO-30 code 



developed by Caughey and Jameson [8] which solves the full-potential equation 
in a conservative, finite-volume form for flow over wing-fuselage combin- 
ations. The wake is treated as a contact discontinuity across which the 
pressure is continuous but the tangential velocity component is discontinuous. 
This inviscid code was also employed by Streett [1,9] to interact with an 
integral boundary-layer method to develop a code called TAWFIVE and further 
details of the procedure are provided in references 8 and 9. 

For the present method, viscous effects are only determined on the wing. 
The inviscid-viscous interaction is performed by running the inviscid code and 
using the inviscid velocity distribution to drive the boundary- layer code. 

The viscous displacement surface calculated by the boundary-layer code is 
added to the original wing surface and the inviscid solution is then deter- 
mined for the new wing surface (wing plus displacement thickness). The dis- 
placement surface updates are slightly under-relaxed to minimize oscillations. 
The number of these "global" iterations varied from about five to twelve 
depending on the Mach number, sweep angle, and angle of attack. 


GOVERNING BOUNDARY-LAYER EQUATIONS 

The nonorthogonal curvilinear coordinate system is illustrated in figures 
1 and 2 and the unit vectors and metric coefficients for this body-oriented 
system are derived in Appendix A. The independent surface variables are s r 
and y-| e where s r refers to s measured along the root-chord line and 
y le refers to y measured along the wing leading edge. The angle 9 
is the angle between the s- and y-coordinate lines. The three-dimensional 
laminar and turbulent boundary-layer equations are derived in detail for this 
coordinate system in Appendix B. For convenience, they are simply listed here 
as: 

continuity equation, 
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spanwise momentum equation. 
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where h^ and h 2 


h 2 are the metric coefficients, Kj and K 2 are the 
qeodesic'curvatures of the curves y le = constant and s r = constant, 
respectively, and K 12 and K 21 are geometric parameters which are functions 
of the metric coefficients and geodesic curvatures. Their explicit values are 
provided in the Appendices. 


The boundary conditions on the above equations are. 
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In order to close the system of equations, some assumptions for the 
Reynolds stresses must be made. For the present study, the Cebeci -Smith Liu J 
eddy-viscosity S model is employed. This relates the turbulent stresses to the 
mean velocity and total enthalpy profiles by. 
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where Pr^. is the turbulent Prandtl number and the 
is defined in the inner layer by [5], 


eddy-viscosity coefficient 



and in the outer layer by [5], 
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... In order to solve eqs. (l)-(4) an inviscid solution has to be DrovidPd 
either by an inviscid code or experiment, and both initial and side^ounda™ 
conditions must be determined since the equations are of the Daraboli! 

coordinate. Therefore, along the wing tip the governing equations reduce to! 
continuity equation. 
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spanwise momentum equation. 
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energy equation. 
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The boundary conditions for the inf inite-swept-wing equations are the same as 
those given for the general equations. 

The root chord line of the wing is assumed to be a line of symmetry where 
the spanwise derivatives th.*S.SriSi 

is d singular^° U However, ^differentiating the ^^ioe momentum 
equation with respect to y 1e yields a nonsingular equation, 
along the root-chord line the governing equations become. 
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where u 5 means the derivative of u with respect to s r . The boundary 
conditions for the attachment- line equations are: 
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TRANSFORMATION OF THE GOVERNING EQUATIONS 
Before proceeding to descritize the governing three-dimensional boundary- 
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Next, define the vectors Z and C by. 


Z = [q] , C = [g] (36) 

Then, using eq. (36), eqs. (27) and (28) may be combined and rewritten as the 
vector equation. 
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Along the wing tip the spanwise derivatives of the flow variables are set 
equal to zero. Thus, the inviscid and geometric terms which are different 
from the general three-dimensional case are: 
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The boundary conditions for 
eq. (40). 


the above are the same as those given by 
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The boundary conditions for the attachment- line equations are the same as 
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2 "2 ^ie **2 ^i7 ^ 5 


- Y-j) ] 


h l Y 2 


m 8 " h l Y 2 2 ls~ ^ K 2 csc0 ^ 

r 
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m 


11 


m 


12 


- 1 + ^7 * h l K 12 Y l + h l Y l 2 35; ,K 2 CSC0) 
= - h L K 2 cote y 1 2 /t 2 


( 57 ) 


°l = w 


C 1 u ref 


°5 = °3 G 


(i - W ) 


(58) 


where the terms y^ - Yg are » 


Y 1 = u 


se 


u 

y 2 = u 


ref 

se 


^1 du se 
3 u se ^ Y le 


Y 4 = u ref W{ e 


h l 5v e 


'l a 


T *- 

Applying the transformation given by e,s. (231-128) to the eddy-viscosity 
formulas given by eq. (7) yields, 

{e + ) = (^i.) /Re f- [0.4 c 2 dr)} 2 [1 - exp (- -J)] 2 


(59) 


m’ • ' v 


*[© 2 + © 2 + 2 h C0SS ^^ ]2 


(60a) 


{e + ) = (^i) /Re (0.0168)|J e c 2 [/I + + 2 ^i cos0 

' m' '•v ; o 


15 



(60b) 


" /F + (P 2 G)^ + 2 7 2 cosgTg ] dr)| 
where Re is a local Reynolds number and. 


A " ^ ^i w c 2 


(^) "eV c 2 do t(S 2 * ft, 2 («)' 


+ 2 h cose (|£) (|6) ]' 

w w 


finite-difference technique 

accurate “<"* • second-order 

scheme is implicit in the normal direction and exp^-ci^H^hr*™ 0 • 6] ’ The 

~- e - - issxzs s,Nri deri ' 

notation - Z ( s r ,y-j e , tj^ ) is employed, where, 


S„ = s + A s 
"HI r i 


i = 1,2,.. . IMAX 


' le j +1 ^le . + Ay i e » J = 1.2....JMAX 


\ + A \ 


k = 1,2, . . .KMAX 


A coordinate stretching procedure develoDed hv Rinti-n 0 ^ rioi • • 

normal direction so that the mesh is c Sstered rnl f ) 15 “ sed ,n the 
Central difference operators are de?iSed by the reU^ons ^ ^ 9 * 


6 h Z 1 > J > k = (Z i,j,k+l ' Z i,j,k-l )/(AT1 k + A Vl } 


y i.j,k ^ Z i,j+l,k " Z i, j-l,k ^ 2Ajr i 
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Applying the above operations at the predictor step to egs. (37), (38), and 
(31) they become, respectively: 


where. 



W/( 



■O 

c 

03 



oo 



(75) 


+ A [a, 

vi L 1 


11 5 i + V 2 » J * k 


A ti G i+ V 2 ,J» k ^ " ° 


The inviscid and geometry parameters are evaluated at the (i+l/jj.j) 

SSTTSiSS*- - SWW" 

equations. 




An 


S.w * 


x-l 


l z j,k * z j,k-l' 


(76) 


- a m Z j,k-1 + b m z j,k • C m Z j,k+1 T u m 'J.k 


+ <L C 4 L = e, 


(77) 


m 


a e E J,k-l + 


b e E j,k 


' c e E j,k+1 b e 


(78) 


where a 






e> 


D e > c e> 


and d p are scalars. 


dp are scaiars, b m and u m 
vector. Their values are provided in 


d m are 2 by 2 


sr s3.“i5 1 *.:’ rv.™ s.s;:;‘.s™ - 

Thomas algorithm [14]. 

The inf inite-sweDt-wing equations, root-chord line equations and 

^fpS^s* Sf'SSffX “«?. ^ ones’ given in ’ 

Appendix C so their values are not listed here. 

The character of the three-dimensional boundary- layer equations was shown 

by Raetz [15] to R a ae P t rfo 0 Jo C rmula U trMrwen-rnown y zo^ro? ’influence and 

5^^’co^ t^Th concept simply stated ^‘ bab , a through 

point p is instantly felt by al P“’ nba P atsinS through that 

p and then convected downstream a ® 9 varies across the boundary 

s ws? ° f * 

solution 1 can ll ldv]nced torn* £'«•*!»• maximum 

step size must be determined from the relation [16], 


As < Minimum {Ay. -A- jr*1 
r all j,k 11 1 


•} 


(79) 


For the present -thod, the spanwise step 

the marching step size is calculated a g proceeds spanwise, along the 

attachment 6 1 i ne, 3 f ronTthe^oot^o the tip using the attachment- line equations. 


19 


are marche<Tstrea^ise°fro^^ three-dimensional equations 

the upper surface trSllJ "ed ,*a d tTe n f ^ihJrtSrhS l^ 9 ^ t0 
lower surface trailing edqe. The sol.itinn I?l th L U hment llr ? e to the 
is obtained in the same manner by solvina thp rnnt G h r °^- dnd tl P' chord lines 
infinite-swept-wing equations resoect i ,2 ^ 0t - cho , rd llne equations and 

conditions for the three-dimensional equations^ Pr ° Vldes the S1de ^ndary 


BOUNDARY-LAYER PARAMETERS 


2 i 

C = s 

f s p u 2 (80a) 

OO OO 

and a spanwise skin-friction coefficient by, 

2 i 

C = y 

p T 7 ( 80b ) 

OO CO 


Using the transformation given bv eas 
ficients become. 


(23)-(28) the skin-friction coef- 



(81a) 


2c, 


^ r Pe u U e u ref -i r &G> 
{ o~ )[ 7~> W 


y /Re p 


W 


(81b) 


The limiting streamline angle, p w , is defined 


as 


P w = tan -1 (C f /C f ) 
y s 

Therefore, the equation for the surface streamli 
given by 


(82) 


nes in the (s,y)-plane is 



= tan p 

w 


(83a) 


20 



so that in the (x.y)-plane the above becomes 


fdyi = r cosQ cosAi t p 
^ L COSK J w 

dx 


(83b) 


w 


A streamwise and a spanwise velocity integral thickness are defined as. 


6 = / (i _ _P H_) dn 

s o p e u e 


(84a) 


pv 


\ = / u p v 
y o e 


-) dn 


(84b) 


In transformed coordinates the velocity integral thicknesses become. 


— [// c ? dT) - 


s /Re " 0 2 


(85a) 


6* = -i- [j" e c 2 dn - ^ 9 (r, )] 
y /Re o e 


(85b) 


The three-dimensional displacement surface. A*, was defined by Moore [17] to 
be the solution of the equation, 


* [P e A * " ^ (p p ** " P ^ dn -' = ° 


e e 


9s. 


★ * -| 

[p u h 0 sin 9 (A - 6 )J 
L e e 2 s 




le 


[p v h. sin 0 (A - 6 )] = 0 
L e el y 


( 86 ) 


When the displacement surface has been calculated for both the upper and lower 
surfaces of the wing, a new wing surface is defined by adding the displacement 
surface to the original wing coordinates, and the inviscid solution is then 
determined for the new wing surface. The displacement surface is not 
smoothed, but is slightly under-relaxed during the first few viscous updates. 
This procedure (global iteration) is repeated until the viscous updates are 
less than some user-specified tolerance. On the plane of symmetry the 
boundary- layer parameters are given by. 


/Re F °° 00 


(87a) 


w 
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(87b) 


C f = 0 

y 


6 s = ~ U 6 c 2 dT1 - ft\)] 
/We o e 


V^[/\ «n-^ 9 (n e )] 

J /Re o v ye e 


dnd along the attachment line they are given by. 


C f =0 
s 


V X/ 2 

C - ° - r e 


P e 'i r u ref u se^ r dG^ 


= 2 Cl “ ^ > ( -T^’ «f> 


6 s = ("IT ”) [/ 6 c 2 dT1 - f{T1 e n 

se o e 


V2 


6 = 


V h 1 ^ h 

C— u ) [/ e c 2 d-n 

se 0 


ref 


9 (n e )] 


(88a) 

(88b) 

(89a) 

(89b) 

(90a) 

(90b) 


DISCUSSION OF RESULTS 

The inviscid solution is determined on a grid which consists of 160 arid 
points in the streamwise direction, 32 points spanwise and 24 points normal to 
the wing-fuselage surface. For the viscous solution a total of 110 streamwise 
and 31 spanwise grid points are used for both the upper and lower surfaces of 
the wing. The number of grid points across the boundary layer varies from 26 
at the attachment line to about 75 at the trailing edge. Each of the 
calculations presented here was performed on a Gould PN 9005 computer which 
has an operating speed of 0.8 MFL0PS. 


As a test case for the present method, a solution was calculated for flow 
over a swept and tapered wing attached to a cylindrical fuselage. The winq 
has a aspect ratio of eight and a NACA 0012 cross-section, without twist. The 
planform is given in figure 4 and has a taper ratio of 0.5 and a leadinq-edge 
sweep angle of 20 degrees. The solution was determined for a Mach number of 
M® * °- 85 W1 th the wing at an angle of attack of two degrees and a Reynolds 
number based on the mean aerodynamic chord of nine million. The results of 
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the present method are compared to those obtained using the TAWFIVE computer 
code developed by Streett and Melson [9], and to the inviscid results from 
FLO- 30. For simplicity, the present method will henceforth be referred to by 
the acronym TIBLT (for Transonic Interactive Boundary-Layer Theory). 

In order to directly compare the two inviscid-viscous interaction 
methods, the boundary layer was assumed to transition from laminar to 
turbulent flow at the first marching step beyond the attachment line. Both 
methods employ 

FLO-30 as the inviscid solver, but the viscous routine in TAWFIVE is an 
integral method while TIBLT solves the boundary- layer equations in the 
differential form. Twelve viscous updates were performed by both interaction 
methods and the number of inviscid iterations was the same for all three 
calculations . 

Figure 5 compares the streamwise distribution of the pressure coefficient 
determined by the three calculation methods. The effect of the boundary layer 
moves the location of the upper surface shock wave towards the leading edge, 
reduces the strength of the shock, and decreases the static pressures in the 
trailing-edge region. The location of the shock wave moves forward as the 
wing tip is approached and near the tip the shock location is nearly identical 
for both of the interaction methods. This may be due to the fact that both 
methods make an inf inite-swept-wing assumption at the wing tip. Due to the 
location of the shock wave, upper surface boundary layer separation was 
predicted at l = 0.505 for the TAWFIVE calculation and at t = 0.510 for TIBLT. 

The determination of boundary layer separation is different for the two 
interaction methods. Since TAWFIVE is an integral boundary- layer method, the 
criterion for separation is when the compressible shape factor reaches a 
critical value of 2.8 [9]. On the other hand, in TIBLT separation is indi- 
cated when the streamwise skin-friction coefficient goes through zero. 

Neither criteria can be said to be the absolute definition of three- 
dimensional boundary- layer separation, however, they may be thought of as 
indicators that separation is imminent [18-21]. When separation is 
encountered in the TIBLT calculation the displacement surface is calculated 
from the separation line to the trailing edge by assuming the terms in the 
brackets in eq. (85) to be constant and equal to their values at separation. 
The lift and drag coefficients computed by each of the calculation methods is 
presented in table 1. Both the lift and drag were slightly larger for the 
TIBLT calculation than for the TAWFIVE calculation. 


Table 1- 

Lift and 

Drag Coefficients 

for NACA 0012 Wing 


FL0-30 

TAWFIVE 

TIBLT 

C L 

0.302 

0.248 

0.275 

C D 

0.0172 

0.0175 

0.0181 
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As a further test for TIBLT, a subcritical and a supercritical Mach 
number solution were calculated for the ONERA M6 wing for which experimental 
results are available [22]. In both cases the angle of attack was 3.06 
degrees and the freestream Reynolds number was 11.7 million, based on the mean 
aerodynamic chord. The M6 wing has an aspect ratio of 3.8, a taper ratio of 
0.56, and a leading edge sweep angle of 30 degrees. The planform is shown in 
figure 6. The subcritical Mach number was M^ = 0.699 and a comparison of 
the streamwise distribution of the pressure coefficient determined by TIBLT 
and TAWFIVE with the experiment of reference 22 is given in figure 7. Viscous 
effects are moderate for this case and both interaction methods compare 
favorably with the measured values. 

Five viscous updates were performed by both interaction methods. The 
computer run time for the TAWFIVE calculation was 5.64 CPU-hours while the run 
time for the TIBLT calculation was 21.87 CPU-hours, approximately 3.9 times 
that of TAWFIVE. This difference may be accounted for by noting that in the 
integral method the shape of the velocity profiles is assumed a priori, while 
in the differential method they are determined as part of the solution. The 
additional calculations performed across the boundary layer at each marching 
step for the differential method results in the longer run times. 

Typical velocity and total enthalpy profiles from the TIBLT calculation 
are provided in figure 8. The spanwise velocity profile exhibits velocity 
crossover; the velocity component is positive near the wing surface and 
negative at the edge of the boundary layer. This type of velocity profile can 
not be modeled by the TAWFIVE code, nor by most integral methods and illus- 
trates an advantage of differential methods over integral methods. The 
locations where this type of spanwise profile occurs on the upper surface of 
the wing is shown in figure 9. The dark triangles in the figure are the 
places where the spanwise profile exhibits velocity crossover. The majority 
of the spanwise profiles along the root chord line have velocity crossover. 

The chordwise extent of this type of profile diminishes as the wing tip is 
approached, and the locations are concentrated near the trailing edge. 

The supercritical solution was for a freestream Mach number of M^ = 0.84 
and a comparison of the streamwise distribution of the pressure coefficient is 
presented in figure 10. Both interaction methods predict a weaker shock wave 
than the experiment and the computed locations of the shock are slightly 
forward of the experimentally observed location. The TIBLT code predicted a 
shock location somewhat closer to the experiment than the TAWFIVE code, except 
near the wing tip where, as in the NACA 0012 solution, the shock location is 
nearly identical for both interaction methods. Surprisingly, the TAWFIVE code 
predicted upper surface boundary- layer separation at l = 0.209, while no 
separation was predicted by the TIBLT code nor was there any indication of 
separation in the experiment. 

Ten viscous updates were performed by both interaction methods and the 
computer run times were 7.82 and 44.53 CPU-hours for the TAWFIVE and TIBLT 
codes, respectively. Figure 11 shows the locations where the spanwise 
velocity profile has velocity crossover for the upper surface of the wing. 
Figure 11 is very similar to figure 9 except that the chordwise extent of the 
crossover profiles is smaller for the higher Mach number case. The angle of 
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attack and the Reynolds number were the same for both cases, thus the 
effect of increasing the Mach number does not change the distribution of the 
crossover profiles but does reduce the number of such profiles. 

The present method was also applied to a Douglas Aircraft Company (DAC) 
wing-fuselage semispan model which was tested by Spaid and Roos [23,24] at the 
NASA-Ames Research Center 14-Foot Transonic Wind Tunnel. The wing has an 
aspect ratio of 6.8, a taper ratio of 0.3 and a quarter-chord sweep angle of 
35 degrees. It employs 11.71 degrees of twist and the planform is shown in 
figure 12. Two experiments were conducted. The first was for a subcritical 
Mach number of M - 0.5, an angle of attack of six degrees, and a Reynolds 

number based on the mean aerodynamic chord of 3.4 million. The second test 

was for a Mach number of M = 0.825, an angle of attack of four degrees, and a 

Reynolds number of 4.5 milTion. A total of 378 static pressure orifices were 
located at nine spanwise stations on the primary wing. An auxiliary wing with 
a two-degree-of-f reedom probe traversing unit was attached to the fuselage 
downstream of the primary wing. The probe traversing unit could be installed 
at any one of the nine spanwise locations on the auxiliary wing, thus allowing 
boundary- 1 ayer profiles to be obtained along each row of static pressure 
orifices. Spaid defined a local skin-friction coefficient and an inviscid 
yaw-plane flow-direction angle as. 


2 x. 


C f = 


p e u Te 


(91) 


P e = tan' 1 [v e /u e ] 


(92) 


Thus, if p is positive, the inviscid flow is moving outboard towards the wing 
tip It was intended to compare the results of both TAWFIVE and TIBLT with 

the experiment, however, the TAWFIVE code computed the displacement thickness 
on the lower surface of the wing to be unrealistically high and the solution 
became unstable. The TIBLT code had no difficulty in computing a solution for 
either Mach number and the results are presented below. 


A comparison of the calculated and experimental upper surface pressure 
coefficient for the subcritical Mach number is given in figure 13. There is a 
very sharp pressure peak near the leading edge at the inner span stations. 
These large inviscid gradients had to be smoothed out slightly d ^ in 9 the 
first few boundary-layer calculations to keep the surface flow attached. A 
total of five viscous updates were performed for this case and required 2 . 
CPU-hours to run. The inviscid flow-direction angle is compared in figure 14 
which shows the inviscid flow is moving inboard towards the wing root over 98% 
of the chord due to the wing twist. Figure 15 is a comparison of the local 
skin-friction coefficient and figure 16 compares the streamwise vel ocity 
inteqral-thickness distributions. The experimental skin-friction coefficient 
was determined by the Clauser chart technique [25] and the integral thickness 
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from the measured streamwise velocity profiles. The agreement between the 
calculated and experimental skin-friction coefficients is excellent at each of 
the span stations. The calculated streamwise velocity integral thickness 
agrees very well with the experiment near the wing tip, however, at the more 
highly loaded mid-semispan stations it is significantly under-predicted. This 
trend was also observed by Spaid [26] when he compared his experimental 
results to various two- and three-dimensional boundary- layer calculations. 
Figure 17 is a reproduction of figure 7 from reference 26 which shows the* 
comparison of the experimental streamwise velocity integral thickness with the 
three-dimensional boundary- layer method of Cebeci, et al. [5]. The solid line 
in figure 17 is Cebeci' s three-dimensional calculation and n is the percentage 
of the semispan. Cebeci used the experimental pressure distribution to drive 
his boundary- layer code but still under-predicted the experimental integral 
thickness at the mid-semispan station. The difference between the experiment 
and the boundary- layer methods is probably due to the transition location 
along the span and flow unsteadiness in the experiment. 

In a recent paper, Cebeci and his co-workers [27] showed that moving the 
transition location has a significant effect on the displacement thickness 
distribution near the trailing edge. For the present method, the flow 
was assumed to be fully turbulent just prior to the location of the upper 
surface pressure peak, which occurs at about £ = 0.015. Cebeci' s calculation 
determined the transition location to be at £ = 0.025. A challenge for 
three-dimensional boundary-layer methods is to properly account for the 
spanwise variation of the attachment line and transition location. At 
present, TIBLT determines an average constant-chord attachment line and 
assumes transition location. In the experiment, a boundary-layer trip was 
located on the upper surface at £ = 0.06, which is aft of the sharp pressure 
peak so that they may have been trying to trip an already turbulent boundary 
layer. The precise location of the experimental transition line and its 
spanwise variation is not known, and may account for some of the observed 
difference between the boundary-layer methods and the experiment. Also, there 
was some flow unsteadiness in the experiment due to the mixing of the tunnel 
flow with the outside air which produced some probe vibrations on the 
auxiliary wing [23,28]. 

The locations where the spanwise velocity profile has velocity crossover 
on the upper surface of the wing is illustrated in figure 18. Once again the 
majority of the spanwise profiles near the root chord are of this type and as 
the wing tip is approached their chordwise extent diminishes and is 
concentrated near the trailing edge. The calculated inviscid and surface 
streamlines for the upper surface of the wing are presented in figure 19. The 
inviscid flow moves inboard as it goes over the wing due to the wing twist. 

The surface flow also moves inboard until just aft of the mid-chord line where 
it then begins to move slightly outboard, except near the root and tip where 
it continues inboard. 

A comparison of the upper surface pressure coefficient for the super- 
critical Mach number is provided in figure 20. The shock location calculated 
by TIBLT is slightly forward of the experimental location at the inner span 
stations. At the outer span station no shock wave is present, probably due to 
the wing twist and to tip effects. A total of ten viscous updates were 
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performed and the calculation took 43 43 .~rs.^ 

inclination angle, skin-friction coeff 21> 22, and 23, respectively, 

thickness distributions are compared in f g simi1ar t0 those from the 

The results for this supercritical case are y velocity profile has 

subcritical case. The u^nP^urface^f the wing is presented in figure 
velocity crossover on the upper sur e of t 9 for the subcritical case, 

24. The number of these prof es 15 a ^ 1ate d inviscid and surface 
surface 5 flow* 4 near ST «hich is moving farther inboard. 


CONCLUSIONS 

A method for determining the loped. an The U effec? t of 

boundary-layer flow over arbitrary sw pt 9 accounted for through a 
the boundary layer on the inv sad ^^ 10 ^ ^ invisci d-viscous interactions 
displacement surface inter ^’° • varies depending on the flow 

conditions? °BateS InThe results of the present study the following 
conclusions are drawn: 


too severe, — . 

(2) The finite-difference t ®^"JJaq e U over t other V f inite-difference 
SesClnatThf ToJyl difference formulas are forced independent of 
fhe Sign of the crossflow velocity component. 

and 

supercritical Mach numbers. 

,4) in each of the cases considered, the spanwise 

exhibited velocity crossover somewhere i^n^t wing planform a re dependent 

this type of velocity profile occurs 9 An a< j va ntage of the 

on Mach number, Reynolds “U i s 9 that this type of spanwise velocity 

Tome TX Sr,t of the solution without any special 
accommodat ion. 
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APPENDIX A. DERIVATION OF THE UNIT VECTORS AND METRIC COEFFICIENTS 


In order to write the governing boundary-layer equations in a body- 
orientated coordinate system, it is necessary to relate the curvilinear 
coordinates to a reference orthogonal coordinate system. This will yield 
expressions for the curvilinear unit vectors and metric coefficients. A 
constant-chord line (figure 1) is defined as. 


£ = 



(Al) 


In general, the wing sweep and cross-sectional shape are functions of £ 
and y . It can be shown that the differential of eq. (Al), holding 
£ = constant, yields. 


tan A = tan A^ e - £A 


(A2) 


where. 


A = tan A. - tan A. 

le te 


(A3) 


The equation for the wing surface is given by. 


F(x,y,z) = z - cf (£,y) = 0 


(A4) 


Therefore, 


,7 = [5L] i ♦ [£ ] j ♦ [SL ] k 

dx ay dz 


dF 


r df ac-i : f 7 „ f af a£ . af n : . r 

L-c sr — J i - Lf — + c u* ~ + — )J j + k 

^ ax ay ay ay 


A A A 

[- tan <] i + [fA + tan A tan k - c -—] j + k 

ay 


(A5) 


and 
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K 


(A6) 


5f 2 1/2 

| VF l = sec < {1 + [tan A sin k + cos < (fA - c - )] 1 


Define an angle r such that, 


tan r = tan A sin k + cos < [fA - c J 


Then a unit vector which is normal to the surface of the wing is defined by 
the equation. 


- _ vf 

e " IvfI 


= -cos 


r sin k i + sin r j + cos r cos < k 


Referring to figure 1, a streamw 
defined by. 


ise unit vector which is orthogonal to e n 


e = cos < i + sin < k 


a 

Next, define a unit vector e which is 
the 7=0 plane (figure 2) by. 


tangent to an t - constant line 


e = sin A i + cos A j 

yo 


Then a spanwise unit vector which is orthogonal 'o JJe region! 

the 7=0 plane in the direction of e yQ may De aer . 


e • e = cos 
y yo 


Thus, 


e y = cos Q e yQ - sin Q k 


= cos Q s 


in A i + cos Q cos A j - sin a k 


Since e was defined to be normal to 
which yields. 


e their dot product must equal 

n* 


( A7 ) 

( A8) 
is 

( A9) 

in 

(A10) 

in 

(All) 

(A12) 

zero 
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f A 1 3 ) 


tan Q = cos A sin r - sj n a cos r sin k 


cos r cos < 


coord i"nate ^system 5 i s 'nonor^hogona ^because U the Tot *"‘d H ° WeVer * the - 
not zero and is given by, use the dot Product of e s and e is 


eos 9 = e $ • e y = cos k cos 0 (sin A - tan < tan a) (A1 

The position vector of a point (TUT) -in • 

,y,z) in Cartesian coordinates is given by, 

^ * i + y j + 7 k A 

f A 1 

(s!y!n)? P However\ sincTT^nf y" ^rv'a^on 1 C00rdinate s ^tem by 
to choose as independent variables (s v ^ n) 9 h hG Wlng ^ ls more conv ^nient 
the root chord and y,„ refers to s f refers t0 s <"ong 

POiht (7,y ,T) may be Uniquely described by (s v le ?n fn9 H* 96 ' ThUS any 
of the position vector on the sorface of the Singes givIX d1ffere "tid' 


d? = ^ dy 'e) ? + d V + «y„) 3 

+ dSr * 5^7 dy le> k ’ 


hq ds e + h„ dv p 
1 r c s 'Z ^le e v 


(A16) 


wing (s,y,0) "to the inLjend^r^rfacrvariabLff ^ ^T) P ° lnt ° n the 
one may obtain, ace variabl es (s r ,y le ,0). From eq. (A16) 


h i -i 

r r r 


h o - • a 9y 7 Qp 

2 y "^37 ^7 J+ ^; k 

Performing the indicated differentiation qi 


gives. 
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(A18a) 


dx 

T5T 


ay 

W. 


dz 

as. 



0 


( — ) cos ic tan k 
v c r 

r 


and 


5x 



= tanA cosQ. cosA, 
le le 


= cosQ, cosA, 
le le 


(A18b) 


(A18c) 


(A19a) 


(A19b) 


dz 


w 


le 


-tanQ cosQ, cosA, /cosA 
le le 


(A19c) 


Thus, the solution of eq. (A17) yields the metric coefficients as, 
c cos K 


u * 

n l ~ C r COS K 


( A20) 


_ cos A le cos Q le 
h 2 COS A COS Q 


(A21) 
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APPENDIX B. DERIVATION OF THE NONORTHOGONAL BOUNDARY-LAYER EQUATIONS 


With expressions for the unit vectors and metric coefficients in the 
curvilinear coordinate system the conservation equations may be derived. The 
velocity vector is given by. 


^ = u e„ + v e + w e„ (Bl) 
s y n 

where u, v, and w are the velocity components streamwise, spanwise and 
normal to the wing surface, respectively. The substantial derivative of the 
velocity vector is therefore. 


D$ _ Du ; 
Ut ~ El s 


, Dv * . Dw 

Uf e y Uf 


De De De n 
e n + u T5T- + v If/ + w UT- 


( B2) 


Thus, in order to calculate the derivative of the velocity vector, one must 
first determine the derivatives of the curvilinear unit vectors. To obtain a 
unit vector along any coordinate line, differentiate the position vector (eq. 
(A15)) with respect to the particular coordinate and divide by its 
magnitude. Thus, 


4r' l|-l 


aft , 

‘ V 


+ (f-) L * 

r 


r az 

l "5s" 



aft , 


(B3a) 


In a similar manner, 

* _ aft , . 
e v = 1 h 2 


(B3b) 


Differentiating eq. (B3a) with respect to y-j e and eq. (B3b) with respect 
to s r yields. 
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a 2 $ 


d 2 ft 


ay le 5s r dy le 


■ Ijrr ( h i e s> = ^7 ^ V = as r oy ie 


(B4) 


Then it follows that, 


&h 2 a ^e y 

wz 


ah. a 5e s 

1 e + h. 


Taking the dot product of each term in eq. (B5) with e $ gives. 


ah ah 9 

ae a i , on i _ 2-i 

y . e = * L— L- - COS 0 -jrr-J 

■Ss; * ^ ^ **r 

A A 

where eq. (A14) was used for the dot product of e s and e y' 
^product of each ter, in eq. (B5) with e y to get. 


( B5) 


(B6) 


Next, take the 


de 
5y le 


. dh 9 ^h 

s • e = i- («=■ - cos 9 






’y n ! ' os r ^le 

Differentiating eq. (A14) with respect to s r yields. 


(B7) 


A 5e w a ae . ae 

Q . y + e . — = - sin 0 tt- 
e s -5C y ^7 ^r 


(B8) 


and using eq. (B6) in eq. (B8) one may obtain, 

* ah. 

•sr- • e v = -R7 th 2 cos 0) ' ^ 


■&s r ‘ y "2 ~~r 


L - 5 s; v 2 

Similarily differentiating eq. (A14) with respect to y, e 9 i'ves, 


ae. 


"s oy 


+ e 


le 


ae s . 00 
• -* — = - sin 0 -gr— 
y ^7 y le 


(B9) 


(BIO) 
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and using eq. (B7) in eq. (BIO) one may obtain, 

A 

5e 

y A 1 rd r dh, 

^7 ' 65 " ^ [ 5yj7 (h i cos ®) - ^ 

Now consider the derivative of p W1 > h 
= s r ,y le ,n) which may be written as* ° S ° me ar 9 ument » Q 


(Bll) 


de_ 


*T = a s e s + b s e y + C S e n 


(B12) 


Taking the dot product of each term 

dot product of each term with e 
which are, y 


in the above with e and then taking the 

will yield expressions for a^ and b 

^ s 


b = cscre (. 


de s 

T5T 


e y } 


de. 


a s = ' cote csce (-5 q 1 * e ) 


(B13) 


Substituting eq. (B13) into eq. (B12) yields 



Following a similar procedure for the derivative of e y it 
ae &e 

y _ f y % r 2 A a 

~5q~ “ ^‘Sq - * V >- csc 0 e - cot 0 esc 0 e 1 + c e 

b y y n 


(B14 ) 

can be shown that, 
(B15 ) 


thPi> e pyniiV ? + Cy ° nly contnbute t0 the normal momentum equation and 
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(B16a) 


ae iA ®b. * o a 

■gs^ = Ti^ ^ h 2 cos0 ^ “ '5y~^' cot0 csc0 e s + csc 0 e y^ 


+ e„ 
s n 


ae 


ah 


dh 


. Ulln UMl *■ 2. \ 

— = i- [■*?=■ - COS0 ^_](-cot0 CSC0 e q + CSC 0 e v ) 


"l 


le 


+ C r e n 

s n 


and using eqs. (B6) and (Bll) in eq. (B15) gives. 


d6 i fth* 5h« o a * 

ai; = ^ [ ^; - cose i^ csc 0 e s - cote csce e y } 


+ c e 
y n 


ae . a 

_JL = jL [_ (h 1 cos0) 

# »le ’T We 1 


ahp A A A 

(csc^0 e $ - cot0 csc0 e^) 


+ c e 
y n 


The substantial derivative for steady flow in the curvilinear coordinate 
system is given by. 


d a ^ a . ,, a 
w =u ¥ v ^y w ¥ 

- u5 + v 5 + w t- 


Substituting eqs. (B16) and (B17) into the above yields. 


(B16b) 


(B17a) 


(B17b) 


( B18) 
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(B19) 


De - de 

W * K 1 “ (- cote e s + csce e y) + + * “ 3F-] e n 


ah 


ah, 


y 1 2 j a a 

+ h"“'h 2 sine ^ ' cos0 ^y^ _cot0 e s + csc0 e y) 


and 


De y U r dh i dh 2i , 

TTt = hj h 2 sine " cos0 "Ss^ ( cscG e s " cot0 e y) 


a a &G A 

+ K 2 v ( csc0 e s - cote e y) + [ c y + + w -5n^] e n 


( B20) 


where Kj and K 2 are the geodesic curvatures of the curves y-| e = constant, 
and s r = constant, respectively, and are given by. 


Sh, 


1 ? \ Ur, 1 

K 1 = h. h ? sin0 hr ^ h 2 cosG ) " 

12 r -Me 


( B21 ) 


dh. 


K 2 = TTp^ 


TTTnS ^ ( h l cos0 ) - 


( B22 ) 


Using eqs. (B19) and (B20) in eq. (B2) gives the substantial derivative of the 
velocity vector as. 


jj|r = - K 1 COt 0 + K 2 CSC 0 + Kj 2 uv] e c 


_ ^2 ^otO v + csc© u + K 2 i uv] e 


+ t TFE + w + W L ' i + (c s + + e n 


de de 


+ w 


De n 

W 


( B23) 


where the parameters and K 21 are given by. 
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1 


( B24) 


12 h^ h 2 sin c 0 UJ 1e 

dh 


[ (l + cos 2 e) -2 cose ^!] 


5h, 


i dn 9 ? 1 1 

K = l— T [^1 (1 ♦ COS 0) - 2 cose ^-] 

di h^ sin 0 r ,e 

Next, consider the gradient of a scalar, p, which may be written as, 

Vp - P s e s + Py e y + P n e n 

The differential of a scalar is given by, 

dp = Vp • dft 

where the differential of the position vector is, 
dft = h x ds r e $ + h 2 dy le e y + dn e n 
Using eqs. (B26) and (B28) in eq. (B27) yields, 

dp = (p $ h 1 + Py h lC os0) ds r + (p s h 2 cos0 + p y h 2 ) dy ]e + P n dn 

However, dp inay also be written as, 

d P = -ST* ds r + WT dy le + M dn 
ua r le 

Equating eqs. (B29) and (B30) and solving for p s and p y gives. 


_ esc 0 2>p _ cote csc0 bp 


p s = 1T—'5s: 


*2 ^le 


esc 0 bp cotQ cscQ 5p 


p y = fi 7 ~ ^y 


le 




( B25) 


(B26) 


( B27 ) 


( B28) 


( B29) 


( B30) 


(B31a) 


( B3 lb ) 
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With eq. (B31), eq. (B26) becomes. 


tfn = r csc 0 d P cot0 csc0 3 d i A 

p r ]e s 


4> r csc 0 dp COt0 CSC0 3 d i * an A 

^ ' — KJ V e y + ^ % 


[14] Ire uIg?«ulI e b r o1I 0 I“ce a 5 "), E " er9y eqUat, '°" S Wr1tte " ,n tensor notatt< 

D\f 3x. . 

p w = " Vp + ^ r. . ’ 1- ’j’ k = i.2,3 ( 

J 

u 2 

ST [ pu j H - FF fe- ( H - 4-) - U 1 T ,-j] - 0 ( 

J J J 

where H is the total enthalpy, p is the dynamic viscosity coefficient 
Pr is the Prandtl number, u T is the magnitude of the total velocity, 

2 2 2 

Uj - u + v + 2 uv cos0 ^ 

and T.j is the stress tensor which is given by, 
r du, 3u • ? 3u. 

T ij = p ^-5xT + Sjr) - T 6 ij ^ (I 


The boundary- layer assumption is that the boundary layer consists of a 
thin viscous layer very close to the surface of the wing [14,29]. Therefore 
the variation of the pressure and the unit vectors across the boundary layer 

^L a th Ummed t °-5 e n ? 9ll91 ’ ble * However, since the velocity components range 
from the inviscid values at the outer edge of the boundary layer down to zero 

tLn h thn ln9 Sur ^ce, the velocity gradients normal to the wing are much larger 
than those parallel to the surface. Thus, in the viscous and heat-conduction 
terms, only the normal derivatives are retained. This produces a parabolic 
c 6 k Rations which may be solved by a streamwise marching procedure 
Substituting eqs. (B23) and (B32) in eq. (B33) yields for the streamwi e 
momentum equation, c 

pu 3u pv 3u A au 2 ? 

Fj" 3s^ + FJ pw m " p K 1 cotG u + P l<2 esc 0 / 
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( B37 ) 


CSC 0 dp COtQ CSC9 ap__ + Q (u gi) 

+ P K 12 uv = ■ ^ ^ M 


and for the spanwise momentum equation, 

pu av + pv &v + pw - p k cote v‘ + pK^ cscO u 


W +, " w * 


av -■ — 2 •• — 2 


+ p k 21 uv "Sy-jg 


csA <*_ t cot^csce |E_ ♦ •_ („ $ 


'1 


( B38) 


Using eq. (B18) in eq. 


(B34) gives the energy equation as. 


pu aH 


+ 


pv dH 


aH 
pw -5n 


a 


[■“ 


W 


aH 

an 


p (i 



a 

W 



( B39) 


Now 

by. 


consider a line element, <U, «hich in the tensor notation [30] is given 

3 3 ( B40) 


dr = I l g n - j dXj 

i=l j=l 

be wr1tte " aS ’ 


dJt 2 = dft • dft = h^ ds r + h^ dy-| e 


+ 2h x h 2 cos 0 ds r dy 1e + dn 


( B41 ) 


Comparison of eqs. (B40) and (B41) gives the elements of 9i j, 


' i J 


"1 

h 1 h 2 cose 


h 1 h 2 cos0 0 


2 


Therefore, g-j j 


"2 

0 0 

i s a symmetric tensor whose determinant is given by. 


( B42) 
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( 843 ) 


Is ( = h x 2 h 2 sin 2 e 


The divergence of a vector in a general coordi 
tensor notation [30] as, y a 1 coord1 


nate system is given in the 


. 1- 1,2,3 

/g t 

where V 1 is the contravariant velocity vector, 

V 1 = V //a — - r u v iT 
V i //g ii - ^ w] 

With eqs. (B44) and (B45) the steady continuity equati 
V • (pfr) = 0 


( B44 ) 


(845) 


on may be derived. 


9 

[/g p V 1 ] = o 


a 

"5s~ 


(puh 2 sinQ] + (pvhjsine) + (pwh^sine) 


(B46) 


Equations (B46), (B37). (B38) anH fR?ai 
boundary-layer equations for a laminar ™ ar « the governing nonorthogonal 
governing tobulent eSuatloL SnL “"Possible flow. To obtain the 

tv ° f - 

uleJ%Je U b^nd t lr?-(^ r 1p^SxlLuor^,^r?1 * 

are multiplied by the normal ve or t ? e f1 “ ctuating quantities which 
because across tL Lnda™ laye? ^ o^ t0 „ 66 r*!" 6 '?- This is 

a-awryni 


pw = pw + p-'w' 


(B47) 


and the governing equations become for a turbulent flow: 
continuity equation, 


gr- (puh ? sin0) 
r c 


^7 


(pvhjSinG) 


d 

5rT 


(pw 


h^sinQ) = 0 


(848) 


42 



streamwise moniGntum Equation, 

r 

pu 9u + pv au + pw 4^- - P K i cote y2 + pK 2 csc9 v 

*" 1 

u au v„ au 2 

re e . e e_ _ K cot0 u c 

+ p K 12 uv = P e L-ry ^ TTJ 1 e 


+ K ? csc0 V 2 + K 12 u e v e ] + In [v- M ' p U w ^ 


(B49) 


spanwise momentum equation, 

p 2 

pu av . pv av + &L . p l cote v + p k 1 csce u 

w 2 


r U e dv e + V e . .cote v 2 
+ P k 21 uv = p e trq- ^ ^ ^ 2 e 

+ K x csce u e 2 + K 2 iU e v e ] + ^ W ‘ P 


(B50) 


energy equation, 

pu SH pv |H_ t p; 9H 


r U- 

tn L Fr 5n 


+ n (1 ' 77^ ^n 


vTTT] 


(B51) 


gradient ten»s_have 


the edge of the boundary layer. 
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APPENDIX C. COEFFICIENTS OF THE TRANSFORMED EQUATIONS 

transformed'variables l°rl 1ayer equations *" 

subscript are evaluated at 1) for £ T iables without an i 

corrector step. The parameter N hafthe va ^ nf S ^ P a ? d + u t li + 1/2) for the 
and two at the corrector step. The terns b h f !5? P r 5 dfctor step 


^m» d m» and e m are defined by. 


b m 


^ml b m3 
b m2 b m4 


, dm = 


d ml d m3 
d m2 m m4 


• e m = 


ml 

; m2 


(Cl) 


The coefficients for the momentum equat 


ions are given by. 


5 j,k- V 2 / 4 Vi 


c " = b j,k+v 2 / n 

b ml ‘ b m4 ■ + S» * > F J,k /As r 



(C2) 


d ml = - m i 0 fF j, k+ i - F j, k _i)/As r 
d m2 = " m i 0 (G j, k+ i ~ G j>|( _ 1 )/As r 
d m3 = d m4 = 0 


t - N 

t l ~ 7 a 2. . 

J,k 
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t - N 

l 2 J 

"7 <9 j + l,k - 9 

j-l,k )/Ay le 



t 3 = 

*10 f 1,j,k /A5 r 




^4a = 

(tj + t 2 + t 3 ) 

* ( F j,k+r 

F j,k-l> 


t 4b 

(tj + t 2 * tj) 

* (G 0,k+l " 

S.k-l 1 


^5a 

*10 (Ar| k + "V- 

■0 F j,k F i, 

,j,k /iS r 


l 5b 

*10 ( fi, V + ln k 

-l> F J.k G i. 

,j,k /4s r 


t 6a = 

- J * 7 K + 

‘Vl 1 G j,k 

( F j*l,k 

- F j-l,k )My le 

II 

_Q 

VO 

+-> 

- J *7 H + 

A Vl ) G j,k 

(G jn,k 

- G j-l,k> My le 

l 7a = 

*2 F j,k + *5 F 

j,k G j,k + 

r 2 

*8 G j,k 

' mn C2 j,k 

t 7b 

*9 F j,k + *4 ' 

j,k G j,k 

.2 

*3 G j,k 

- *12 c 2 j>k 

e ml = 

1 l 4a + C 5a + 1 

6a - 7 K 

+ fi Vi' 

1 t 7a 

e m2 : 

■ l 4b + "Sb + 1 

6b - 7 (*\ 

+ A \-l 

) t 7b 


At the corrector step the e m vector has the additional terms 


(C4) 


(C5) 
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e ml " e ml 


e m2 = e m2 


+ C m (F f,j,k + l - F j,j, k ) - (F fJjk - 
+ C m (G i,j,k +1 - G f,j,k) - \ ( G i,j,k - S.j.k-l) 


(C6) 


For the solution of the energy equation the term f. mea ns f- , =, t 

the predictor step and f, +1 • at the corrector rh! 1+ l/2. j,k at 

coefficients are defined by * 3 ’ corrector Step. The energy equation 


" e " CTl J,k- V 2 ' AVl 


c e = °1, 


J,k+l/2 


/ An, 


b e ' a e + c e + m 10 t 4 \ + 4 Vl> F j,k /4s r 


(C7 ) 


t _ N 
t 8 " 7 °2 . , 

J,k 


= ™10 (f j,k - f i,j,k) /4s r 
l 10 = T m 7 (9j +1 ,k - 9 j-l,k )/4y le 
*11 = * *9 + *10) * ( E j,k t . - E j,k.l) 

* 12 ■ -t m 7 G j,k K + 4 Vi)( E jti,k - E j-l,k) / 4 n e 


tl3 ' N ° 4 j.k-V2 (F l> k - 1 ' F j. k)M ” k -l 


t 14 = N 0 


4 j,k+l/ 2 tFj ’ ktl " F j. k,/iT> k 
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OA unknown point 
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Flgure 3. Finite-difference cell. 
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